/**********************************************************************
First look at Uganda and then at Kenya
********************************************************************************/
set more off
use "${data}\UgandaKenya20162019Replication.dta", clear


********************************************************************************


global did1 1.postban postxnearbordc  // nearbordc -- 0/1 dummy (1 = less than mean)
global did2 1.postban postxlnbordc  // postxlnbordc - continuous variable
global didkm c.dist_bordc##c.dist_bordc##i.postban 
global did3 1.postban postxproximity   // postxproximity - continuous variable
global didkenya 1.postban postxprotect  
global logitcontrols pa dist_close dist_nairo maize_prod ///
		ruggedness pop_2015 dist_dirtr elevation slope_perc ///
		i.year
global twcontrols year yrmaize

eststo clear

replace dist_bordc = dist_bordc/100
replace bordc_sq = dist_bordc^2

replace postxbordc = dist_bordc*postban
replace postxbordcsq = bordc_sq*postban
replace year = year-2016
gen yrmaize = year*maize_prod



*-----------UGANDA-------------------------------------------*

if 0 {
preserve
keep if country == "Uganda"
xtset id year

foreach num of numlist 10 30 50 {
xtreg defor`num' 1.postban 1.postban#1.near_bordc  $twcontrols if gte10`num' == 1, fe cluster(code2)
est sto c1`num'
eststo c1`num', add(RMSE e(rmse))
sum defor`num' if yr < 2018 & gte10`num' == 1
eststo c1`num', add(Mean r(mean))

xtreg defor`num' 1.postban 1.postban#c.ln_dist_bordc  $twcontrols if gte10`num' == 1, fe cluster(code2)
est sto c2`num'
eststo c2`num', add(RMSE e(rmse))
sum defor`num' if yr < 2018 & gte10`num' == 1
eststo c2`num', add(Mean r(mean))

xtreg defor`num' 1.postban 1.postban#c.proximity  $twcontrols if gte10`num' == 1, fe cluster(code2)
est sto c3`num'
tempname dydx


postfile `dydx' marg sedydx zvalues using "${data}\proximityuganda`num'.dta", replace 
                 
foreach z of numlist 2(1)18   { 
lincom   1.postban#c.proximity*`z'/1000 

post `dydx' (r(estimate)) (r(se)) (`z')

}
postclose `dydx'

eststo c3`num', add(RMSE e(rmse))
sum defor`num' if yr < 2018 & gte10`num' == 1
eststo c3`num', add(Mean r(mean))
}
esttab c1* using  "${tables}\UgandaDID.tex", replace ///
se nonotes  style(tex)  b(%12.3f) se(%12.3f) scalars(N Mean RMSE)    sfmt( %12.0gc %6.2f %6.2f) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" ""  ) ///
	keep (1.postban 1.postban#1.near_bordc ) nonumbers fragment  ///
	coeflabel(1.postban "2018-2019"  1.postban#1.near_bordc "2018-2019 x $<$  100 km to border" ) ///
	 prehead({ \begin{tabular}{lccc}    ///
	 \hline \hline ///
	 & \multicolumn{3}{c}{Forest threshold cutoff} \\ ///
	 & 10\% & 30\% & 50\% \\ ///
	& (1) & (2) & (3)  \\ ) ///
	prefoot( \\ ) postfoot( )
	  

esttab c2* using  "${tables}\UgandaDID.tex", append ///
se nonotes  style(tex)  b(%12.3f) se(%12.3f) scalars(N Mean RMSE)    sfmt( %12.0gc %6.2f %6.2f) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" ""  ) ///
	keep (1.postban 1.postban#c.ln_dist_bordc) nonumbers fragment  ///
	coeflabel(1.postban "2018-2019"  1.postban#c.ln_dist_bordc "2018-2019 x ln(km to border)" ) ///
	 prehead( \hline &  (4) & (5) & (6)  \\ ) ///
	prefoot( \\ ) postfoot( )
	  
	  
esttab c3* using  "${tables}\UgandaDID.tex", append ///
se nonotes  style(tex)  b(%12.3f) se(%12.3f) scalars(N Mean RMSE)    sfmt( %12.0gc %6.2f %6.2f) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" ""  ) ///
	keep (1.postban 1.postban#c.proximity) nonumbers fragment  ///
	coeflabel(1.postban "2018-2019"  1.postban#c.proximity "2018-2019 x proximity"  ) ///
	prehead( \hline &  (7) & (8) & (9) \\ ) ///
	prefoot( \\ Grid FE & yes & yes & yes  \\ ///
	Time trend & yes & yes & yes & yes \\ ) ///
	  postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{Each panel represents estimations using a different measure of distance from the ///
	  nearest official Kenyan border crossing. The estimator is a linear probability model with standard errors robust to clustering  at the county-level. All ///
	  regressions include grid cell fixed effects, a time trend, and the time trend interacted with maize ///
	  productivity. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )
	  
	  
/*
*************** CHECK ON OMITTED VARIABLES *****************	  
xtreg defor30 1.postban 1.postban#c.proximity  $twcontrols if gte1030 == 1, fe cluster(code2)

local r = e(r2)*1.3
psacalc delta 1.postban#c.proximity , rmax(`r')
psacalc delta 1.postban , rmax(`r')	  
estimates clear
restore
*/
}

*-----------UGANDA VARIATION IN BINARY IMPACTS DEFINTIONS -----------*

if 0 {
	preserve
	keep if country == "Uganda"
	xtset id year
	replace dist_bordc = dist_bordc*100
foreach num of numlist 10 30 50 {
	forval d=  25(25)400 {
	gen near`d' = dist_bordc<`d'
	gen treat`d' = postban*near`d'
	loc f = 0
	xtreg defor`num' 1.postban treat`d' $twcontrols if gte10`num' == 1, fe cluster(code2)	
		if _b[treat`d']/_se[treat`d'] > 2  {
				loc f = e(F)
				eststo c`num'
				eststo c`num', add(Km `d')
				eststo c`num', add(Fst `f')
			}
		drop near`d' treat`d'
		}
	
	}	

esttab c* using  "${tables}\Ugandathreshold.tex", replace ///
se nonotes  style(tex)  b(%12.3f) se(%12.3f) scalars(N Km )    sfmt( %12.0gc %6.2f %6.2f) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" "" )  ///
	keep(1.postban treat*) nonumbers fragment  ///
	coeflabel(1.postban "2018-2019" treat150 "$<$ `d' km x 2018-2019")  ///
	 prehead({ \begin{tabular}{lccc}    ///
	 \hline \hline ///
	 & \multicolumn{3}{c}{Forest threshold} \\ ///
	 & 10 & 30 & 50 \\ ///
	& (1) & (2) & (3)  \\ ) ///
	prefoot( Grid FE & yes & yes & yes  \\ ///
	Time trend & yes & yes & yes & yes  \\ ) ///
	  postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{Standard errors are robust to clustering ///
	  at the county-level. All regressions include grid cell fixed effects, a time trend, and the time trend interacted with maize /// 
	  productivity. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )
	estimates clear
	restore
	}



*-----------UGANDA SCOBITS-------------------------------------------*
if 0 {
	preserve
keep if country == "Uganda"
xtset id year

*create means for CRE models*
loc mx1 " "
loc mx2 " "
loc mx3 " "
loc mx4 " "

foreach x in  postxnearbordc yrmaize {
	bys id: egen m`x'=mean(`x')
	loc mx1 = "`mx1' m`x'"
}
drop  myrmaize
foreach x in  postxlnbordc  yrmaize {
	bys id: egen m`x'=mean(`x')
	loc mx2 = "`mx2' m`x'"
}
drop myrmaize
foreach x in  postxproximity yrmaize  {
	bys id: egen m`x'=mean(`x')
	loc mx4 = "`mx4' m`x'"
}


foreach num of numlist 30 50 {
	
loc logL = -999999999
forval a=  0.2(0.05)0.7 {
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	scobit defor`num'  i.postban##i.near_bordc  $twcontrols  `mx1' if gte10`num' == 1, vce(cl code2) constraint(1) 
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				eststo c1`num': margins, dydx(postban) at(near_bordc= (0 1)) predict(pr) post noestimcheck
				eststo c1`num', add(alpha `a')
				eststo c1`num', add(logL `ll')
			}
		}
		
	}	
	
	
loc logL = -999999999
forval a=  0.2(0.05)0.95 {
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	cap scobit defor`num'  i.postban##c.ln_dist_bordc $twcontrols  `mx2' if gte10`num' == 1, vce(cl code2) constraint(1) 
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				eststo  c2`num': margins, dydx(postban) at(ln_dist_bordc = (0 1 3 5)) predict(pr) post noestimcheck
				eststo  c2`num', add(alpha `a')
				eststo  c2`num', add(logL `ll')
			}
		}
		
	}	
	loc logL = -999999999
forval a=  0.2(0.05)0.95 {
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	cap scobit defor`num' i.postban##c.proximity $twcontrols  `mx4' if gte10`num' == 1, vce(cl code2) constraint(1) 
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				eststo c3`num': margins, dydx(postban) at(proximity= (0 .25 .75 1)) predict(pr) post noestimcheck
				eststo c3`num', add(alpha `a')
				eststo c3`num', add(logL `ll')
			}
		}
		
	}	

}

esttab c* using  "${tables}\UgandaScobit.tex", replace ///
se nonotes  style(tex)  b(%12.3f) se(%12.3f) scalars(N alpha ll )    sfmt( %12.0gc %6.2f %6.2f) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" "" "" "" "" ) ///
	 nonumbers fragment   ///
		 prehead({ \begin{tabular}{lcccccc}    ///
	 \hline \hline ///
	 & \multicolumn{3}{c}{30\% forest threshold} & \multicolumn{3}{c}{50\% forest threshold} \\ ///
	& (1) & (2) & (3) & (4) & (5) & (6)  \\ ) ///
	prefoot( Grid FE & yes & yes & yes & yes & yes & yes \\ ///
	Time trend & yes & yes & yes & yes & yes & yes  \\ ) ///
	  postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{Standard errors are robust to clustering ///
	  at the county-level. Estimator is a scobit, with grid cell means of all included variables to accommodate fixed effects. The alpha ///
	  for the scobit is defined by a grid search that chooses the alpha with the highest log likelihood function.  ///
	  All regressions also include  a time trend,and the time trend interacted with maize productivity. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )

estimates clear

restore

}
*-----------KENYA ESTIMATIONS-------------------------------------------*

if 0 {
preserve
keep if country == "Kenya"
la var year Year
*create means for CRE models*


foreach x in postban postxprotect yrmaize {
	bys id: egen m`x'=mean(`x')
}


xtset id year
xtreg defor10  1.postban year yrmaize if gte1010 == 1,  fe cluster (code2) 
est sto c1
eststo c1, add(RMSE e(rmse))
sum defor10 if yr < 2018 & gte1010 == 1
eststo c1, add(Mean r(mean))

xtreg defor10  1.postban##1.pa year yrmaize if gte1010 == 1,  fe cluster (code2) 
est sto c2
eststo c2, add(RMSE e(rmse))
sum defor10 if yr < 2018 & gte1010 == 1
eststo c2, add(Mean r(mean))


xtreg defor30  1.postban year yrmaize if gte1030 == 1,  fe cluster (code2) 
est sto c3
eststo c3, add(RMSE e(rmse))
sum defor30 if yr < 2018 & gte1030 == 1
eststo c3, add(Mean r(mean))

xtreg defor30   1.postban##1.pa  year yrmaize if gte1030 == 1,  fe cluster (code2) 
eststo c4
eststo c4, add(RMSE e(rmse))
sum defor30 if yr < 2018 & gte1030 == 1
eststo c4, add(Mean r(mean))


xtreg defor50  1.postban  year yrmaize if gte1050 == 1,  fe cluster (code2) 
eststo c5
eststo c5, add(RMSE e(rmse))
sum defor50 if yr < 2018 & gte1050 == 1
eststo c5, add(Mean r(mean))


xtreg defor50  1.postban##1.pa year yrmaize if gte1050 == 1,  fe cluster (code2) 
eststo c6
eststo c6, add(RMSE e(rmse))
sum defor50 if yr < 2018 & gte1050 == 1
eststo c6, add(Mean r(mean))

esttab c1 c2 c3 c4 c5 c6 using  "${tables}\KenyaDID.tex", replace ///
se nonotes  style(tex)  b(%12.3f) se(%12.3f)  scalars(N Mean RMSE)   sfmt(%7.0fc %4.2fc %4.2fc %7.0fc %4.2fc ) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" ""  "" "" "" "" ) ///
	  keep(1.postban 1.postban#1.pa year)  order(1.postban 1.postban#1.pa year) coeflabel(1.postban "2018-2019" 1.postban#1.pa "2018-2019 x PA" year "Trend") nonumbers fragment  ///
	 prehead({ \begin{tabular}{lccccccc}    ///
	 \hline \hline ///
	 &\multicolumn{2}{c}{10\%}  &\multicolumn{2}{c}{30\%}  &\multicolumn{2}{c}{50\%}   \\ ///  ///
	& (1) & (2) & (3) & (4) & (5) & (6) \\ ) ///
	prefoot( Grid FE & yes & yes & yes & yes & yes  \\ ///
	Time trend & yes & yes & yes & yes & yes & yes \\ ) ///
	  postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{The dependent variable is an indicator equal to one of deforestation occurred in a grid cell in a given year.  The percentages in the header indicate the samples used for estimation, where the percent is the amount of baseline forest cover required to consider the grid cell to have had forest in the baseline.  The estimator is a linear probability model with standard errors are robust to clustering at the ///
	  county-level.  All regressions include grid cell fixed effects, a time trend, and the time trend interacted with maize ///
	  productivity. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )
restore
}

if 0 {
preserve
estimates clear
keep if country == "Kenya"
la var year Year
loc x0 "year yrmaize"
loc x11 "postxproximity yrmaize"
 
global treat i.postban##i.pa

*create means for CRE models*


foreach x in postban postxprotect yrmaize {
	bys id: egen m`x'=mean(`x')
}
	  
	 **** .2 is optimal alpha, but the loop seems not to be working correctly  so re writing this code to just run for .2
loc a = .2
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	cap scobit defor30  i.postban year yrmaize   myrmaize if gte1010 == 1, vce(cl code2) constraint(1) 
	eststo c1: margins, dydx(postban year)  noestimcheck predict(pr) post 
	eststo c1, add(alpha `a')
	eststo c1, add(logL  e(ll))		
sum defor10 if yr < 2018 & gte1010 == 1
eststo c1, add(Mean r(mean))

	cap scobit defor30  i.postban year yrmaize   myrmaize if gte1030 == 1, vce(cl code2) constraint(1) 
	eststo c2: margins, dydx(postban year) noestimcheck predict(pr) post 
	eststo c2, add(alpha `a')	
sum defor30 if yr < 2018 & gte1030 == 1
eststo c2, add(Mean r(mean))
	cap scobit defor50  i.postban year yrmaize  myrmaize if gte1050 == 1, vce(cl code2) constraint(1) 
	eststo c3: margins, dydx(postban year) noestimcheck predict(pr) post 
	eststo c3, add(alpha `a')
	eststo c3, add(logL e(ll))	
sum defor50 if yr < 2018 & gte1050 == 1
eststo c3, add(Mean r(mean))
	cap scobit defor10  i.postban##i.pa year yrmaize mpostxprotect myrmaize if gte1010 == 1, vce(cl code2) constraint(1) 
	eststo c4: margins, dydx(postban year) at(pa = (0 1)) noestimcheck predict(pr) post 
	eststo c4, add(alpha `a')
	eststo c4, add(logL e(ll))				
sum defor10 if yr < 2018 & gte1010 == 1
eststo c4, add(Mean r(mean))
	cap scobit defor30  i.postban##i.pa year yrmaize  mpostxprotect myrmaize if gte1030 == 1, vce(cl code2) constraint(1) 
	eststo c5: margins, dydx(postban year) at(pa = (0 1)) noestimcheck predict(pr) post 
	eststo c5, add(alpha `a')
	eststo c5, add(logL e(ll))		
sum defor30 if yr < 2018 & gte1030 == 1
eststo c5, add(Mean r(mean))
	cap scobit defor50  i.postban##i.pa year yrmaize mpostxprotect myrmaize if gte1050 == 1, vce(cl code2) constraint(1) 
	eststo c6: margins, dydx(postban year) at(pa = (0 1)) noestimcheck predict(pr) post 
	eststo c6, add(alpha `a')
	eststo c6, add(logL e(ll))


esttab c1 c2 c3 c4 c5 c6 using  "${tables}\KenyaScobit.tex", replace ///
se nonotes  style(tex)  b(%12.3f) se(%12.3f)  scalars(N Mean logL  alpha)   sfmt(%7.0fc %4.2fc %4.2fc %4.2fc ) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" ""  "" "" "" "" ) ///
	  nonumbers fragment ///
	 prehead({ \begin{tabular}{lccccccc}    ///
	 \hline \hline ///
	 &\multicolumn{6}{c}{Forest threshold sample}  \\ /// 
	 & 10\% & 30\% & 50\%  & 10\% & 30\% & 50\% \\ ///
	& (1) & (2) & (3) & (4) & (5) & (6) \\ ) ///
	prefoot( Grid FE & yes & yes & yes & yes & yes  \\ ///
	Time trend & yes & yes & yes & yes & yes & yes \\ ) ///
	  postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{Standard errors are robust to clustering at the ///
	  county-level. Marginal effects shown. Estimator is the scobit, with the choice of alpha ///
	  defined by a grid search that chooses the alpha with the highest log likelihood function.  ///
	  All regressions include grid cell fixed effects, a time trend, and the time trend interacted with maize productivity. ///
	   Scobit estimations include variables demeaned at the ///
	  grid cell level rather than fixed effects (the Mundlak adjustment). * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )
	  
restore

}

if 1 {

*------------------------Figures -------------------------------------

use "${data}\proximityuganda30.dta", clear

gen LB95 = marg-invttail(100000,.025)*sedydx
gen UB95 = marg+invttail(100000,.025)*sedydx
 
replace zvalues = 1/(zvalues/1000)
la var marg "Marginal effect"
la var LB95 "95% CI"
 
twoway connected marg LB* UB* zvalues, xtitle(Km to border) ytitle(Increase in probability of deforestation)  legend(order(1 2  )) ///
lpattern(l - - ) lwidth(medium thin thin)  msymbol(i i i) lcolor(cranberry gs10 gs10)
graph export "${figures}\deforuganda.eps", replace
}